Excised black hole spacetimes: quasi-local horizon formalism applied to the Kerr 
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We present a numerical work aiming at the computation of excised initial data for black hole 
spacetimes in full general relativity, using Dirac gauge in the context of a constrained formalism 
for the Einstein equations. Introducing the isolated horizon formalism for black hole excision, we 
especially solve the conformal metric part of the equations, and assess the boundary condition 
problem for it. In the stationary single black hole case, we present and justify a no-boundary 
treatment on the black hole horizon. We compare the data obtained with the well-known analytic 
Kerr solution in Kerr-Schild coordinates, and assess the widely used conformally flat approximation 
for simulating axisymmetric black hole spacetimes. Our method shows good concordance on physical 
and geometrical issues, with the particular application of the isolated horizon multipolar analysis 
to confirm that the solution obtained is indeed the Kerr spacetime. Finally, we discuss a previous 
suggestion in the literature for the boundary conditions for the conformal geometry on the horizon. 
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I. INTRODUCTION 

Trying to accurately describe black holes solutions as 
evolving physical objects in numerical simulations is of 
direct interest in astrophysics. Numerical simulations 
have made a great leap forward in the past few years, 
mainly with the first stable simulations of black hole 
mergers in full general relativity by Pretorius [1], Cam- 
panelli et al. [2] , Baker et al. [3] , and a few other groups 
(see [4] for a review) . Several of these simulations model 
black holes in their equations by punctures. These punc- 
tures basically change the topology of spacetime to han- 
dle evolution of singular objects (see [5] for the first pro- 
posal of this method). 

Notable exceptions to this are references [1, 6-8], 
that use an excision approach for black hole evolution: 
the simulations only evolve the spacetime outside two- 
spheres that are supposed to encircle the black hole sin- 
gularities. 

In the same spirit, we are here trying to describe black 
holes as physical objects characterized by their horizons. 
Defining the physical laws for event horizons of black 
holes has been notably done by [9], [10] and [11] in a 
so-called membrane paradigm. The aim was to describe 
them as fluid-like two-membranes with physical proper- 
ties. However, applying evolution laws to event horizons 
is problematic due to their teleological (non causal) be- 
havior (see for example [12] for a description of the phe- 
nomenon). Being defined as a global property of space- 
time, a local notion of causality actually docs not apply 
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to event horizons. 

Alternative local characterizations have been formu- 
lated in the past 15 years by [13], [14] and [15]. They are 
based on the concept of trapped surfaces, dating back to 
Penrose's singularity theorem [16]. Defined locally, those 
objects behave in a causal way in a general dynamical 
context, with local evolution laws following from the pro- 
jection of Einstein equations, e.g the Navier-Stokes [17] 
"fluid bubble" analogy or the area evolution law in [18]. 
For this work we shall use local characterizations for iso- 
lated horizons, prescribing the physics of non-evolving 
black hole horizons. 

Following the prescriptions of [19], [20], [21] and pursu- 
ing the numerical explorations of [22, 23] and [24] among 
others, we try to numerically implement those objects as 
boundary conditions imposed on the (3+1) form of Ein- 
stein equations, in a three-slice excised by a two-surface 
(single black hole case). This is done here using the fully 
constrained formalism (FCF) of [25] (see also [26]), with 
maximal slicing and Dirac gauge, based on the (3+1) 
formalism, and with spectral methods-based numerical 
resolution using the LORENE library [27]. An impor- 
tant point here is that we drop out the usual conformal 
flatness hypothesis and solve for the conformal geome- 
try, so that we can exactly recover a slice of a stationary 
rotating vacuum spacetime. 

Contrary to free evolution schemes, which are the most 
used prescription for (3+1) simulations in numerical rel- 
ativity, an important feature of constrained schemes is 
the necessity to solve constraints on each three-slice, in 
the form of elliptic equations. These equations generally 
need additional conditions to be imposed on grid bound- 
aries, following reasonable geometrical and physical pre- 
scriptions. Our approach particularly requires a specific 
handling of boundary conditions for the two dynamical 
gravitational degrees of freedom. This is a crucial point 
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of our calculation; we will justify and apply here a no- 
boundary treatment for these quantities. 

The paper is organized as follows: we first review 
in Sec. II fundamental geometrical properties associated 
with isolated horizons in general relativity. In Sec. Ill, 
we quickly give the basic features of the fully constrained 
formalism and the methods we use to treat the confor- 
mal part. Sec. IV discusses implementation of boundary 
conditions for the system of equations, and specifically 
discusses the conformal metric part. Sec. V gives the 
numerical results obtained, and confronts them to a bat- 
tery of tests characterizing the physics of the solutions. A 
discussion follows in Sec. VI, in regard of previous works 
concerning the computation of the conformal part in the 
black hole initial data problem. We also raise the ques- 
tion of applicability of this scheme to other more general 
astrophysical cases. 

Throughout all this paper, Greek letters will denote 
indices spanning from to 3, Latin indices from k to m 
shall denote indices from {1,2,3}, and indices from a to 
c have the range of {2,3}. All formulae and values are 
given in geometrical units (G = c = 1). We also use the 
Einstein summation convention. 



II. ISOLATED HORIZONS AS A LOCAL 
DESCRIPTION OF BLACK HOLE REGIONS 

A. Trapped surfaces and expansion 

The concept of a trapped surface in a Lorentzian frame 
has been first defined by Penrose in 1965 [16] in connec- 
tion with the singularity theorems. It relics on the no- 
tion of expansion of the light rays emitted from a surface, 
that we explain here. We start by a closed spacelike two- 
surface S embedded in spacctimc, topologically related 
to a two-sphere. We assign to it a two-metric q a b in- 
duced by the ambient four- metric g^, and its associated 
area form ef 6 . Two future null directions orthogonal to 
S are associated with this surface. Representative vector 
fields are denoted I th and being respectively oriented 
outwards and inwards. We can for example assume that 
our spacetime is asymptotically flat, so the orientation 
can be defined without ambiguity. 

The expansion 9g of S along £ M is the area rate of 
change along this vector: Ci ef fc = 9g ef 6 . £ is the Lie 
derivative, here along the vector i^. Same definition goes 
for the vector field k^. For a two-sphere embedded in 
a Minkowski spacetime (flat metric), we have typically 
Ok < and 9i > 0. S is said to be a trapped sur- 
face if both expansions are negative or zero: 9g < , 
Ok < 0. This clearly characterizes strong local curvature. 
A marginally (outer) trapped surface will be character- 
ized by 9 1 = and Ok < 0. Two theorems make the 
connection between those objects and black holes: pro- 
vided the weak energy condition holds, the singularity 
theorem of Penrose [16] ensures that a spacetime con- 
taining a trapped surface necessarily contains a singu- 



larity in its future. Following this result, provided the 
cosmic censorship holds, another result by Hawking and 
Ellis [28] conveys that a spacetime containing a trapped 
surface necessarily contains a black hole region enclosing 
this surface. 

Marginally outer trapped surfaces (MOTS) are in- 
tended as models for the black hole boundary (see [29] 
for a discussion of its relation with the boundary of the 
black hole trapped region). Local horizons (trapping 
horizons for [13], isolated and dynamical horizons in [14], 
[15]) are defined as three-dimensional tubes "sliced" by 
MOTS, with additional geometrical properties. The iso- 
lated horizon case is detailed below; for a review, sec [30] . 

In vacuum stationary spacetimes, all horizons are at 
the same location, which is also the location of the 
event and apparent horizons (constructed with outermost 
MOTSs). In the more general case, and assuming cosmic 
censorship, local horizons are always situated inside the 
event horizon in general relativity. 



B. Isolated horizons 

The notion of isolated horizon is aimed at describing 
stationary black holes. It is based on the notion of non- 
expanding horizons, and defined as a three-dimensional 
tube H foliated by MOTS, and with a null vector field V 1 
as generator. The three- metric induced on the tube has 
then a signature (0, +, +). 

We also define a (3+1) spatial slicing for our spacetime, 
and S a two-slice of our isolated horizon at a certain value 
of the time parameter t. The spacelike two-metric on S 
is denoted q a b- 

The shear tensor a a b on the two-surface is defined along 

as 

Cab = ~ [£( qab - 9 e q ab } ■ (1) 

Using the fact that 9g = and the dominant energy 
condition, the Raychaudhuri equation for null tubes [21] 
ensures that both the shear along V 1 and the energy- 
momentum tensor projected on i^, and evaluated on the 
surface must vanish: a a b = and T a b£ a £ b = 0. These are 
additional properties constraining the geometry of the 
horizon. 

An isolated horizon is also required to be such that the 
extrinsic geometry of the tube is not evolving along the 
null generators: [Cm,Dj] = 0, where is the con- 
nection on the tube induced from the ambient spacetime 
connection. If this last condition is dropped, we only re- 
trieve a non-expanding horizon (linked to the notion of 
"perfect horizon" [31]). 

The isolated horizon formalism has already been stud- 
ied extensively in numerics, as a diagnosis for simulations 
involving black holes, where marginally trapped surfaces 
are found a posteriori with numerical tools called ap- 
parent horizon finders (See for example [32], [33], [34], 
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[35]). Let us note that very often, apparent horizon find- 
ers actually locate MOTS on the three-slice considered 
(not necessarily outermost ones). A thorough study of 
geometrical properties of isolated horizons located a pos- 
teriori can be found in [36]. 

In the present paper we employ isolated horizons as 
an a priori ingredient in the numerical construction of 
Cauchy initial data for black hole spacetimes. More pre- 
cisely, we impose conditions on the excised surface char- 
acterizing it as the slice of a non-expanding horizon (see 
below). This approach to the modeling a black hole hori- 
zon in instantaneous equilibrium has been investigated in 
[19-24, 38], where a prescription for the conformal met- 
ric is assumed. The main feature of this work is the 
inclusion of the conformal metric in the discussion, not 
through analytical prescriptions, but indeed by numeri- 
cal calculation. This problem has also been recently ad- 
dressed in [39]. In Sec. IV we tackle the description of our 
special treatment of the conformal metric on the excised 
boundary and compare with previous results, in particu- 
lar through the numerical recovery of excised Kerr initial 
data. 



III. A FULLY CONSTRAINED FORMALISM 
FOR EINSTEIN EQUATIONS 

All the following is a summary of the physical and 
technical assumptions set in [25]. For a review of (3+1) 
formalism in numerical relativity, the reader is referred 
to [40], or more recent reviews like [41] and [42]. 



stein equations system(see for example [43]): 



R + K - KijK 13 



DjKl 



D,K = 



(4) 
(5) 



dt 



K, 



C B K., 



-DiDjN 



N {R l3 - 2K ik K$ + KKij} ,(6) 



Di and Rij being, respectively, the connection and the 
Ricci tensor associated with the three-metric 7^ . Quan- 
tities without indices represent tensorial traces. These 
equations are referred to respectively as the Hamiltonian 
constraint, momentum constraint and evolution equa- 
tions. 



B. Conformal decomposition, maximal slicing and 
Dirac gauge 

Now we must choose a set of variables and a gauge, to 
get a partial differential equations system that we solve 
numerically. The first ingredient in the formalism pre- 
sented in [25] is the conformal decomposition of the three- 
metric [44]. We define on each slice an extra metric 
noted fij , that will have a vanishing Riemann tensor (flat 
metric) and will be time independent. The existence of 
such a metric in a neighborhood of spatial infinity is en- 
sured by our sub-manifold being asymptotically flat. The 
associated flat connection is noted T>i. We introduce in 
E t a conformal metric such that its determinant coincides 



with that of fij, as: 



A. Notations and (3+1) decomposition 



(7) 



We consider an asymptotically flat, globally hyperbolic 
four-dimensional manifold Ai, associated with a metric 
9iiv of Lorentzian signature (—,+,+,+). We define on 
Ai a slicing by spacelike hypersurfaces E t , labeled by a 
timclike scalar field t; in this way, the four-metric can be 
written in its usual (3+1) form: 

g^dafdx" = -N 2 dt 2 + llj {dx l + (3 l dt){dx j + f3 j dt), (2) 



The tensor field h lJ we use to encode the conformal de- 
grees of freedom is the deviation of the conformal metric 
from the flat one: 



(8) 



We also define in our equation sources the following 
conformal traceless extrinsic curvature: 



here N and fi l are the usual lapse scalar field and shift 
vector field. 7 y is the spacclikc three-metric induced on 
E t . 

We also define the second fundamental form of E t , or 
extrinsic curvature tensor, as: 



3 



(9) 



We choose for a gauge the generalized Dirac gauge for 
the conformal metric: 



1 



(3) 



with the future-directed vector field normal to E t . 
Writing the vacuum Einstein equations with this formal- 
ism, one comes up with the classical (3+1) vacuum Ein- 



V k j kl = V k h Kl = 



ki 



(10) 



and we add to this prescription the maximal slicing con- 
dition, i.e the vanishing of the trace in the extrinsic curva- 
ture: K = 0. Therefore, A lJ contains all the information 
about extrinsic geometry. 
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Under those conditions, we can rewrite the (3+1) Ein- 
stein equations in what we shall call the FCF system: 

AV; = S^(N^,(3\A^,0), (11) 
A(Ni') = S m) {N^,^,A»,0), (12) 

A/3' + h&VjF = Sl(N, i>, 0), (13) 

,•)-/,'.• N 2 ()],■■■> 

-£r-W* h " - 2 ^ + W« = 

S^(N,i,,^,AV,h"). (14) 

A is the usual scalar flat laplacian (which expression 
from a spectral point of view is recalled in the Ap- 
pendix A). The actual sources S..., in general non-linear 
in the variables and time-dependent, can be retrieved by 
the reader from [25]. 

We must supplement this system with the kinematical 
relation between the three-metric and extrinsic curvature 
of the slice, deduced from (3) and (9) (see equation (92) 
of [25]). This fully constrained scheme is strictly equiva- 
lent to the one presented in [25]. A slightly different ver- 
sion has been presented recently in [26] , focusing on non- 
uniqueness issues. Although the scheme in [26] would 
probably pose no additional difficulty in the present set- 
ting (except maybe some more boundary conditions to 
prescribe to additional variables), there has been no sig- 
nificant indication of problems involving non-uniqueness 
of solutions in our study, that suggested modifications of 
the original formalism. 

Here we choose as variables the quantities Nip, ip, (3 l 
and h lJ . We especially come up with three elliptic equa- 
tions, two scalar and one vectorial. Those are derived di- 
rectly from the Hamiltonian and momentum constraints 
of the (3+1) system, together with the trace part of the 
dynamical equations. In an evolution scheme, these will 
be the conditions enforced at each value of time t. We 
do this for one particular slice. 

We are then left in general with a second-order ten- 
sorial hyperbolic equation (14) dealing with the variable 
/i u , that is obtained by the geometrical relation between 
7ij and Kij, and the dynamical part of Einstein equa- 
tions. 

The goal here is to simulate as accurately as possible 
stationary spacetimes containing one black hole, repre- 
sented by an isolated horizon. In this respect, we shall 
assume a coordinate system that is adapted to stationar- 
ity This will mean that a stationary timclikc Killing vec- 
tor field will be identified with our time evolution vector 
field (Jj)' '• Using this prescription, all the time deriva- 
tives in our equations vanish, so that our sources and 
operators simplify somewhat. In particular, the tensorial 
equation is written using an ellipticlike operator acting 
on : 

Ahl3 ~ ^2^ hlJ = S 2 , N, /?, A ij ). (15) 

Our problem is then totally equivalent to an actual initial 
data problem, where quantities have to be determined on 



a three-slice by elliptic equations, before evolving them. 
The main difference with classical initial data schemes 
like the conformal transverse traceless (CTT), the ex- 
tended conformal thin sandwich (XCTS) scheme or the 
conformal flat curvature (CFC) system, is an additional 
elliptic equation for the conformal geometry of the three- 
slice. Up to now, a vast majority of initial data compu- 
tations have been done using an ad hoc prescription for 
the conformal geometry. The most common one is the 
conformally flat approach, where ylj is simply approxi- 
mated to be the 3D flat metric. This has been done in 
numerous computations, and this type of initial data is 
the most frequently used for black hole evolution simu- 
lations. However, though this conformally flat approxi- 
mation turns out to be well-behaved in most cases, we 
know that it is a strong limitation when trying to com- 
pute stationary black hole spacetimes: it has been proven 
that the rotating Kerr-Newman spacetime does not ad- 
mit any conformally flat slice (see [45, 46]). 

Other prescriptions for the conformal geometry include 
data suggested by the post-Newtonian formalism [47], 
or superposition of additional gravitational wave content 
(see [48]). Let us mention the work of [49] for neutron- 
star binary initial data, which also computes the confor- 
mal geometry using a prescription in [50] , that considers 
as well the dynamical Einstein equations for the confor- 
mal variables. Finally, the exact scheme we have expli- 
cated above has been applied by one of the authors in the 
case of a single rotating neutron star in equilibrium [51]. 
It has led to the computation of strictly stationary initial 
data, that can be directly extended into future and past 
time directions. This is exactly what we are trying to do 
here in the black hole case. 



C. Resolution of conformal metric part 

Apart from the boundary condition problem (that we 
discuss in Sec. IV), our approach for the resolution of 
the tensorial equation presents some peculiarities that 
we explain here. 

The system of equations is composed of equation (15) 
and the gauge condition: 

T>ih ij = 0, (16) 

that we supplement with a condition on the determinant 
of 7 IJ , following from our definition of the conformal fac- 
tor: 

det(f j ) = det(h ij + f j ) = 1. (17) 

We are left with a tensorial equation for a symmetric 
tensor with four constraints: the system has two degrees 
of freedom. We now try to make them explicit and solve 
for the related variables. Any second-rank symmetric 
tensor h lJ can be decomposed in the following way into 
a divergence free part, and a symmetrized gradient part: 

h ij = V l W j + T> j W i + h% , (18) 
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with T>ih 1 ^ = 0. We shall use here variables associated 
only with the divergence- free part , meaning that the 
gauge component (gradient part) of the tensor considered 
has no influence on them. We choose to encode the infor- 
mation in hj, in the two scalar spectral potentials A and 
B presented in [52] , and whose definitions arc quickly re- 
called in the Appendix A. (A more extensive study shall 
be performed in [53]). 

What is remarkable about quantities A and B is that 
they can actually be decomposed into scalar spheri- 
cal harmonics, and that the tensorial Poisson equation 
A/i y = decouples into scalar elliptic equations A and 
B (see the Appendix A) . This is not exactly the case for 
the nonlinear modified elliptic operator (16); however, in 
our numerical scheme, we just slightly modify the sources 
of the equation at each iteration so that we can write: 

AA-^CpCpA = A s (h ij ,N,^,l3,A^), (19) 
KB-^CpCpB = B s (h^,N,^,p,A^), (20) 

the elliptic operator A being defined in the Appendix A. 
We keep the Lie derivative notation for scalar fields, 
to show that this is directly related to the operator in 
Eq. (15); of course, in the scalar case, this operator sim- 
ply reduces to CpA = (3 l T>iA. During the iteration of the 
resolution algorithm, the sources of the equations are up- 
dated so that they stay coherent with the original equa- 
tion in h l i . Equations (19) and (20) are the two elliptic 
equations that we solve at each iteration. 

Specifically, at each step, we proceed as follows: once 
the scalars A and B are determined by the resolution of 
(19) and (20), the Dirac gauge and unit determinant con- 
ditions allow us to totally reconstruct a divergence-free 
tensor, as the expected solution of our tensorial equation. 
This is done by inverting two differential systems ((B5) 
and (B6) of the Appendix B), that express the Dirac 
gauge conditions and definitions of the scalars A and B 
in function of the tensor components. Those differential 
systems involve scalars, which are components of h 13 in 
a tensor spherical harmonics basis (see the Appendix A, 
and Sec. V of [52]). 

The differential systems require three boundary con- 
ditions on the excised surface, to be inverted (see [53]); 
we discuss them in Sec. IV, in a detailed description of 
the scheme. In addition, the trace of our tensor with 
respect to the flat metric is kept fixed, so that the cal- 
culated determinant at this step is one. This reduces to 
an algebraic nonlinear condition for the tensor compo- 
nents. Finally, we update our sources for the next step. 
We note here that the resolution for the variable C intro- 
duced in the appendix is not necessary in this scheme: for 
a divergence-free tensor, C is unambiguously determined 
by the knowledge of B and the trace. 

With this tensorial scheme, the gauge is necessarily 
enforced by construction, so no gauge-violating mode 
can occur. This is in the same spirit as the global fully 



constrained formalism (equations (11-13)) for our equa- 
tion system, that forbids a priori all constraint-violating 
modes. We also emphasize the fact that, in the gen- 
eral case and with an arbitrary source for (15), we do 
not recover an actual solution of the equation by recon- 
structing our tensor this way. This is only true if the 
elliptic equation admits a solution that actually satisfies 
the Dirac gauge and the determinant condition. We can 
see it as an integrability condition for our equation, that 
is for example not generically true during an iteration. 
However, since in our case we are looking for stationary 
axisymmetric data for a single black hole, we know that 
our entire system does admit a solution: it is the Kerr- 
Ncwman spacetime in Dirac gauge. As a consequence, if 
our scheme converges, we know that the tensor field h lJ 
we obtain shall satisfy the dynamical Einstein equations, 
thus equation (16). 

The missing ingredient for solving all our system of 
equations in an excised spacetime is the setting of bound- 
ary conditions for our partial differential equations, fol- 
lowing part of the geometrical prescriptions of the iso- 
lated horizon formalism, namely non-expanding horizon 
boundary conditions. 



IV. BOUNDARY CONDITIONS AND 
RESOLUTION OF THE FCF SYSTEM 

A. Boundary conditions for the constraint 
equations 

Besides the prescription of asymptotic flatness at infin- 
ity and the bulk stationarity prescription, all the physics 
of our system will be contained in the boundary condi- 
tions we shall put on our excised surface. This section 
follows largely the prescriptions of [21]. 

We consider our excised two-surface to be a slice of 
an non-expanding horizon, i.e. a MOTS with vanishing 
outgoing shear. Following Sec. II, this translates into 
several geometrical prescriptions, namely the vanishing of 
the outgoing expansion and the shear two-tensor: 0i = 
and a a b = 0. 

Being an instantaneous non-expanding horizon, the 
evolution of the excision surface will be a null tube. Since 
we are adapting our coordinates to stationarity, another 
important condition on the excised boundary consists in 
prescribing the time evolution vector field of our coor- 
dinates to be tangent to the null tube. Thus, we are 
ensured that our horizon location stays instantaneously 
fixed during an evolution. Those prescriptions on the 
horizon will suffice to give four boundary conditions for 
the constraint equations (one is scalar and the other vec- 
torial), as we see below. 

We certainly have freedom to prescribe the coordinate 
location of our excision surface in our coordinate system. 
For simplicity, we choose the surface to be a coordinate 
sphere, fixed at a radius t-h- We shall denote by s l the 
unit outer spacelikc normal to the surface, that will be 
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tangent to the three-slice E t . The shift vector is then 
decomposed into two orthogonal parts adapted to the 
geometry of the horizon: (3 l = bs l — V 1 . 

The vanishing of the expansion can be expressed as a 
condition for the conformal factor on the horizon: 



4? Di ln(V>) + D l s l + V" 



0, 



(21) 



where wc have used the conformal rescaling s l = i(j 2 s 1 , 
and the notation Di for the connection associated with 
the conformal three-metric. Multiplying (21) by ip, it can 
be seen as a non-linear Robin condition for the quantity 
ip. The requirement for the time evolution vector field 
on the horizon to be tangent to the null tube provides 
the equality b — N. This is a natural way to fix the 
component of the shift normal to the two-sphere. Wc 
must also fix V 1 , the part of the shift tangent to the two- 
surface. For this we make use of the vanishing of the 
symmetric shear tensor <j a b- It can be shown ([22] and 
[21]) that the vanishing of the shear is equivalent to the 
following equation for V 1 : 



<lh, 



l D a V c + q ac 2 D b V c - q ab 2 D C V C = 0. 



(22) 



Here D is the connection associated with q a i, on the 
surface. This means that V 1 is a conformal Killing sym- 
metry for the two-sphere (in particular, quantities in (22) 
can be substituted by tilded conformal ones). Defining 
coordinates (9, ip) on our two-sphere, we prescribe V 1 as: 



V 



dip 



(23) 



and we shall verify a posteriori that this is a (confor- 
mal) axial symmetry. The constant Q will be called the 
rotation rate of the horizon, (p being the azimuthal co- 
ordinate. In the case of the Kerr spacetimc, there is an 
analytical relation between the arcal radius of the appar- 
ent horizon, the (reduced) angular momentum parameter 
and Q. From a more general point of view, different 
values for f2 will likely affect directly the angular momen- 
tum. In the general case, we define a parameter a for the 
angular momentum associated with the entire spacetimc, 
from the dimcnsionlcss relation: 



J 



K 



M AD 



M 



M 2 

1V1 ADM 



(24) 



with Madm the ADM mass of the 3-slice, and Jk the 
Komar angular momentum of the 3-slice at infinity; the 
latter is tentatively defined with the (presumably) Killing 

vector (see Equations (7.14) and (7.104) in [41] for 

explicit expressions for Madm and Jk)- Note that we do 
not impose any Killing symmetry, except on the horizon: 
we know however, by the black hole rigidity theorem [28], 
that an accurate resolution of Einstein equations would 
impose this vector to be so. We discuss the dependence 
between all those quantities in Sec. V. 



Once we have set boundary conditions for the confor- 
mal factor and the three components of the shift vector, 
we must still fix the lapse function on the horizon. Differ- 
ent prescriptions have been considered in the literature 
(e.g. [19-24]). In the spirit of the effective approach in 
[22, 23], we arbitrarily impose the value of the lapse to 
be a constant N-n on the excised sphere. 

As mentioned before, previous boundary conditions de- 
fine with no ambiguity our excised surface to be a slice 
of a non-expanding horizon. Moreover, our choice for 
the lapse, the horizon location and the conformal Killing 
symmetry on the horizon fixes coordinates on the two- 
surface. Only the conformal two-geometry of the excised 
sphere remains to be fixed. This is done in relation with 
the resolution scheme for the h lJ equation. 



B. Boundary conditions for the h rj equation 

We recall that, with the approach developed in Sec. Ill, 
the resolution of our tensorial problem in Dirac gauge re- 
duces to two elliptic-like scalar equations, to be solved 
on a three-slice excised by a two-sphere: we should nor- 
mally provide two additional boundary conditions for 
those equations. 

A result by [54] shows that in the full evolution case for 
this tensorial equation (equation (14)) and in a Dirac-like 
gauge, the characteristics of the equation are not entering 
the resolution domain when the spacetimc is excised by 
a null or spacclikc marginally trapped tube. This means 
that in the evolution case, once the initial data are set, 
there is no boundary condition whatsoever to prescribe 
to the hyperbolic equation. 

The problem is of course different here, where we are 
left with an elliptic equation instead of a hyperbolic one. 
However, a simple analysis will hint that in our particular 
single horizon case, there will not be any inner boundary 
condition to be prescribed on our data. 

Let us examine the case of the elliptic equation in A, 
that we recall here: 



AA 



CpCpA^A s {h l \N^ 7 P,A^). 



(25) 



We will try and exhibit a simplified linear operator act- 
ing on the variable A, that will contain the most relevant 
terms. The double Lie derivative operator acting on A 
can be separated in: 



V> 4 

j^CpCpA 







Z-ipydjiA+jpiCpCfiAy; (26) 



N 2 



the second term contains all the remaining components 
of the double Lie derivative. 

At this point, and with a fixed system of spherical co- 
ordinates, we are allowed to make a decomposition into 
spherical harmonics for all the scalar variables. We write 
in this respect: 



A = ^Ai m Yi m (6,ip), 

(l,m) 



(27) 
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where Yi m are the spherical harmonics of order (/,m), 
defined as eigenfunctions of the angular Laplace operator: 
A 0lfi Y lm = -1(1 + l)Y lm . 

We now point out the fact that, due to our coordinate 
choice, we have (/3 r )(; = o) = (^?)(J=o) 011 the horizon. 
We can use a second-order Taylor expansion to write the 
I = part of the factor in front of the first term of (26), 
close to our surface coordinate radius r^: 



— (0 r ) 2 
N 2 [P 1 



d 2 A = [l + a(r-r H ) 



(1=0) 

S(r - r H ) 2 + 0(r - r H ) 3 ]d 2 A, (28) 

where a and 6 are two real numbers that can be directly 
computed during one iteration, from the values of N, ip 
and r at the excised surface. Our global equation can 
be rewritten for each spherical harmonic I as: 

d 2 2d 
[-a(r - r n ) - S{r - r H ) 2 } —A lm + -—A v 



dr 2 



dr 



- l -^P^ A im = ^s + ^(C £ A)Z, (29) 

where we only keep on the left-hand side the terms given 
in (28), and put the rest (denoted with **) with the 
source. The latter contains the remaining components 
of the double Lie derivative, and involves either terms 
that are not second-order in the radial derivative, or that 
are multiplied by the higher harmonics of j^(0 r ) 2 (sup- 
posedly smaller than the main term, explicitly developed 
in (29)). Thus, we have isolated a linear operator Q a s, 
depending on two real numbers a and 6: 



QaS 



[—a(r - rn) - S(r 
2 d 1(1 + 1) 



d 2 



r dr 



-I. 



(30) 



other contributions being taken as source terms. This 
operator is different from the ordinary Laplace opera- 
tor by the factor in front of the second order differential 
term, which vanishes on the excision boundary. It can 
be shown that the space of analytic solutions on R 3 mi- 
nus the excised horizon belonging to the kernel of Q a s 
is generally of dimension one. This is in contrast with 
the case of the Laplace equation, where it is of dimen- 
sion two. In practice, this will mean that for a numerical 
resolution of an equation Q a sf = Sf, there is only one 
boundary condition to fix for the unknown, mainly the 
behavior at infinity. No additional information is needed 
at the excised boundary for the effective operators (30). 
Operators of this kind arc known in the mathematical lit- 
erature as elliptic operators with weak singularities [55]. 

The very same scheme can be applied to the equa- 
tion (20) for B, the only difference being that the original 
Laplace operator is replaced by a slightly modified one 
(see the Appendix A). 

As explained in Sec. Ill C, after solving for the two main 
equations (29) and its equivalent in B, the inversion of 



the gauge differential systems (B5-B6) (cxplicited in the 
appendix B) for the reconstruction of h tJ requires three 
extra boundary conditions, in addition to the vanishing 
of all quantities at infinity. We obtain them as com- 
patibility conditions based on the original elliptic tenso- 
rial equation (15): we express three decoupled elliptic 
scalar equations for three components of h lJ in the spin- 
weighted tensor spherical harmonics basis, denoted h rr , r\ 
and //, which are directly related to the usual tensorial 
components of h lJ and defined in Appendix A. From the 
tensor equation (15), we deduce: 



2dr] i 2rj 2h rr V 4 

"! — q — ' 
r or 

Ah rr - 



SI 3 ) (31) 



N 2 



{CpC^Y = (S$y '(32) 



6h r 



4 1h ib A 



ST\(33) 



where (fi,r),rr) superscripts indicate the corresponding 
components of h lJ in the tensor spherical harmonics basis 
(see Appendix A). As for the equation involving A, we 
can rewrite the above equations by extracting the weakly 
singular operator Q a g acting on the principal variable, 
the other contributions being put on the right-hand side 
of the equations. For example, the equation in \i can be 
rewritten the following way: 



2 6> M 2(i 

QaSW + + -o = 

r Or r z 



(34) 

with the Lie derivative term containing all the left-hand- 
side contributions of equation (31) not taken into ac- 
count. We do not need to invert this equation: how- 
ever, as the leading order term in Q a g vanishes at the 
horizon (m = l)j we can accordingly write a Robin- like 
boundary condition for the fx quantity: 



/' 



N~ 2 



Adfi (A 6v + 2) 
r dr r 2 

(35) 

This will be used as a boundary condition for the 
gauge system (B5), the source terms being computed 
with quantities from the previous iteration. Using the 
same method, we can write very similar expressions (that 
we do not explicitly give here) for the fields h rr and 
77, to be used as Robin boundary conditions applied to 
the gauge differential system (B6). The three boundary 
conditions are sufficient to invert the two gauge systems 
(B5-B6) [53], and reconstruct the whole tensor from 
the tensor spherical harmonics components (see the Ap- 
pendix B and [53] for details). 

To summarize, the method employed for the resolution 
of the whole /i y system is iterative and can be decom- 
posed for each step in the following way (more technical 
details are provided in the Appendices): 

1. After calculating the source S2 from equation (15), 
we deduce the right hand side of the equation (29) 
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for A, using values from the previous iteration. The 
same is done for the quantity B and its correspond- 
ing source terms. 

2. We invert equations (29) and its equivalent for B, 
only by imposing that the fields are vanishing at 
infinity. 

3. We compute the value of the trace from equa- 
tion (17) (a more explicit expression can be found 
in equation (169) of [25]). This allows us to 
write the two differential systems (Dirac gauge sys- 
tems) mentioned in Sec. Ill C and expressed in Ap- 
pendix B, involving the spherical harmonics com- 
ponents of (scalar quantities). 

4. We invert these two gauge differential systems using 
three boundary conditions similar to Eq. (35), for 
the three scalar spherical harmonics components 
h rr , r\ and fi. As those are compatibility condi- 
tions expressing information already contained in 
Eq. (15), we provide in this way no additional phys- 
ical information. This gives us the spherical har- 
monics components of h lJ . 

5. We reconstruct the whole tensor from the spher- 
ical harmonics components. 

We have not proven here that no boundary condition 
has to be put generically for the resolution of the two 
scalar equations involving A and B in the tensorial sys- 
tem. However, if we implement numerically the resolu- 
tion by the inversion of the operator Q a g at each itera- 
tion, we will not have to impose any boundary condition, 
but only informations coming from the Einstein equa- 
tions. Moreover, a convergence of the entire h lJ system 
would support the coherence of the reasoning, and hint 
that there is, in our case, a deeper physical motive pre- 
venting the prescription of additional information on the 
horizon. The results in Sec. V show this is the case. 



V. NUMERICAL RESULTS AND TESTS 
A. Setting of the algorithm 

All the numerical and mathematical tools we use here 
are available in the open numerical relativity library 
LORENE [27]. Our simulation is made on a 3D spheri- 
cal grid, using spherical harmonics decomposition for the 
angular part and multidomain tau spectral methods (sec 
[56] for a review). The mapping consists in four shells 
and an outer compactificd domain, so that infinity is 
part of our grid and we have no outer boundary con- 
dition to put at a finite radius. Our grid size is typically 
N r x Ng x N v — 33 x 17 x 1. We also have checked 
our code by setting N v = 4, to verify that no deviation 
from axisymmetry occurred. Our innermost shell has a 
boundary at the radius r-n , which will be the imposed lo- 
cation of a MOTS, and will be used as the unit of length 



in all the results presented here. We impose the values 
of all the fields to be equivalent at infinity to those of a 
flat three-space. Finally, trying to get stationary data, 
we prescribe our coordinates to be adapted to this sta- 
tionarity, so that all the time derivatives in the Einstein 
(3+1) system are set to zero. However, even if we expect 
to get axisymmetric data (the only vacuum stationary 
solution for a black hole being the Kerr solution), we are 
always able to solve our equations in three dimensions. 

We proceed with our scheme in the following way: dur- 
ing one iteration, all the variables are updated immedi- 
ately after they have been calculated, so that the sources 
for the next equations are modified. The tensorial equa- 
tion for is the last solved in a particular iteration, and 
we obtain at each step a local convergence for the whole 
tensorial system (including the determinant condition), 
before we proceed to the update of all quantities, and to 
the next iteration. 

We impose on the sphere of radius r~n the conditions 
of zero expansion (21) and shear (22), via respectively a 
Robin condition on the quantity ip and a Dirichlet condi- 
tion on the partial shift V 1 . We also impose the horizon- 
tracking coordinate condition on the radial shift compo- 
nent b. Having set the shape and the location of the 
surface in our coordinates, we are only left with two free 
parameters, which are the boundary value of the lapse 
function and the rotation rate O. As we said, the lapse 
function, which is merely a slicing gauge choice, is fixed 
to a constant value < Nu < 1 on the horizon. We gen- 
erate two sets of data on our three-slice, spanning the ro- 
tation rate from zero (Schwarzschild solution) to a value 
of about 0.22, where our code no longer converges. One 
set will give the solution for the whole differential system 
(the non-conformally flat (NCF) data, supposed to con- 
verge to the rotating Kerr solution), while the other will 
compute conformally flat (CF) data, by putting = 0. 
From a spacetime point of view, the CF data can also 
be seen as a computation of black hole spacetime using 
the so-called Iscnberg- Wilson-Mathews approximation to 
general relativity [57], [58]. 



B. Numerical features of the code 

Figure 1 presents, on the one hand, the absolute ac- 
curacy obtained for the Einstein constraints (in the form 
expressed in [40]) in the NCF case. Regarding fulfill- 
ment of the Einstein dynamical equation, Figure 1 also 
shows the accuracy of the NCF fully stationary solution, 
as well as its violation in the conformally flat case. We 
see the expected improvement for precision of resolution 
of dynamical equations in the full NCF case. Let us note 
that a verification of the gauge conditions is not even 
necessary, as it is fulfilled by construction (we only solve 
for variables satisfying the gauge). This is one of the 
strengths of our algorithm. 

A non-trivial issue of our computation is the link be- 
tween the two physical characteristics of the system (the 
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FIG. 1: Accuracy for Einstein equations resolution (in the 
original (3+1) version of [40]) as a function of dimensionless 
parameter MuQ (see section VC for a definition of Mn)- 
Data are absolute maximum error values for both constraint 
equations, and dynamical equations in both cases. Data are 
taken with N r — 33, N$ = 17, N v = 1. Lapse on the horizon 
is N n = 0.55. 
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FIG. 2: Dependence on the parameter M«f2 of the ADM 
mass, the Komar angular momentum at infinity Jk and the 
angular momentum parameter -j^ (both defined in section 
IV A) for both cases. The value of the lapse on the horizon is 
here fixed at Nn = 0.55. 



mass and angular momentum of the data) and the two 
input quantities supposed to fix them, namely the bound- 
ary value for the lapse and the rotation rate O on the 
horizon. We choose here in our sequence to fix the value 
of the horizon coordinate radius, removing it from the 
list of variables. Results are shown in figure 2. The value 
of the lapse being also fixed, we observe that an increase 
in Q not only affects the angular momentum, but also 
the ADM mass of the spacetime. Moreover, fixing the 
rotation rate does not amount to the prescription of the 
angular momentum to an a priori given value. A de- 
crease in the value of N-h on the horizon results also in 
an increase in jj (defined in section IV A). This stems 
from the fact that our choice for the slicing directly influ- 
ences in this approach the physical parameters (e.g the 
areal radius) of the solution obtained. We note also that 
for a fixed value of Nn, the correspondence between f2 
and -fj is slightly different in the conformally flat case 
and in the NCF case. With our algorithm, a larger value 
of the lapse gives a slightly better convergence of the code 
for high rotation rates of the black hole (until Nu = 0.8 
approximatively). For each lapse the code stops converg- 
ing at a certain value of the rotation rate. We do not yet 
know whether this is a problem of our algorithm to be 
improved, or if this has deeper physical reasons: constant 
values for the lapse and the rotation rate might not be 
"good" variables for the Kerr black hole in Dirac gauge, 
once we reach high rotation rates. The only conclusion 
we can draw from this is that there is a non trivial cor- 
respondence between our "effective parameters" Nn. and 
O, and the physical ones, namely the ADM mass and Ko- 
mar angular momentum. This correspondence is likely to 
be one to one for values of -jj below a certain threshold 
of about 0.85. Reaching higher values for jj is left to 
future numerical investigations. 



Let us mention again the remark made by [22] about 
the boundary condition for the lapse in the XCTS 
scheme. Although it is necessary to fix the slicing of 
the spacetime by an arbitrary boundary condition on the 
lapse, we have the freedom to decide what kind of con- 
dition to impose. The authors in [22] suggest that an 
arbitrary condition of Neumann or Robin type would be 
preferable, because it is more flexible in view of a numer- 
ical algorithm. In particular, not fixing a value for the 
lapse on the horizon, but rather giving a first order pre- 
scription, allows the data to "adapt" to potentially high 
tidal distortions. However, having also tried to impose 
Neumann conditions for the lapse in our configurations, 
we do not see any clear improvement in the robustness 
of the algorithm. This is why we still keep a Dirichlct 
boundary condition as the simplest prescription. 



C. Physical and geometrical tests for stationarity 

One of the tests of stationarity to be made can be the 
comparison between the ADM mass and the Komar mass 
at infinity, defined with the (presumably) Killing vector 

( J^) 4 (equation (7.91) of [41]). The results of this test are 
displayed in figure 3. The comparison between the ADM 
mass and the Komar mass is actually directly linked to 
the Virial theorem of general relativity put forth by [59] . 
The concordance between those masses is equivalent to 
the vanishing of the Virial integral, and has been also 
used as a stationarity marker by [60]. 

We have also computed in both NCF and CF cases 
an estimate of the amount of gravitational radiation con- 
tained outside the black hole in the 3-slice. Following the 
prescriptions of [14], we calculate the difference between 
the ADM mass and what could be called the isolated 
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FIG. 3: Different diagnostics for stationarity in both cases, 
comparing physical quantities at the horizon and at infinity. 
The viriai error computes the difference between ADM and 
Komar masses at infinity, rescaled with the ADM mass, and 
using asymptotic behaviors of the lapse and conformal factor. 
The radiation energy content outside the black hole (resp. 
outer angular momentum content) is the absolute difference 
between the horizon mass Mn (resp. angular momentum on 
the horizon Jn) and the ADM mass Madm (resp. Komar 
angular momentum at infinity Jk), rescaled with the ADM 
mass (resp. Komar angular momentum at infinity). 



horizon mass, defined in geometrical units by: 



M u = 



^R n + 4 J n 
2R H ' 



(36) 



where Rn is the areal radius of Tt. Mn is nothing but 
the formula for the Christodoulou mass [61] calculated 
from the Komar angular momentum Jn on the horizon 
(defined with the same supposed Killing symmetry as 
Jk)- If we have an isolated horizon extending to future 
infinity, the difference between Madm and Mn gives ex- 
actly the radiation energy emitted at future null infinity 
for the data [14]. In non-stationary cases (for example 
binary systems), this is an appropriate estimate of the 
radiation content at an initial given time. 

Results for comparison between the two cases studied 
here are shown in figure 3. Although the gravitational 
energy available for NCF spacetimes is contained under 
10 -7 whatever the rotation rate might be, in the CF case, 
the increase in energy with jj is patent. This measure of 
energy available with respect to jj gives us a way of ap- 
proximating a priori the amount of what is usually called 
"junk" gravitational radiation, that could be emitted on 
a spacetime evolution with conformally flat initial data. 

In the same spirit, we have also computed the accu- 
racy in the verification of a Pcnrosc-likc inequality for 
axisymmetric data, that can be written as: 

e A = — ^ 4^ < 1, (37) 



MM 2 ADM 



V M ADM ~ J ^ 



and Madm is the ADM mass at infinity. Being a little 
more stringent that the actual Penrose inequality, it has 
been first proposed by [28] for axisymmetric spacetimes. 
This inequality is supposed to be verified for all axisym- 
metric data containing an apparent horizon, and to be 
an equality only for actual Kerr data (this is referred in 
[63] as Dain's rigidity conjecture [62]). The results arc 
presented in figure 4. We observe that, if the equality is 
very well verified in the actual Kerr case, this is definitely 
not true for CF data, even for reasonable values of jj. 
In [63] (cf. [64] for a general context), it has been pro- 
posed that this quantity €a (Dain's number) should be 
understood as a strong diagnosis tool for distinguishing 
between Kerr horizons and other isolated or dynamical 
horizons. This numerical observation shows strong sup- 
port in favor of this claim, pulling apart actual Kerr data 
and reasonable approximations of these data. Let us also 
point out the virtual costlessness of this tool, as we only 
have to rely on a single real value. 

We also note that, when computing the rescaled differ- 
ence of Komar angular momentum between the horizon 
and infinity JK J~ I ^ n , we come up in all cases with a dif- 
ference at the level of numerical precision for resolution 
(see figure 3). This is of course coherent with the fact 
that gravitational waves cannot carry any angular mo- 
mentum in axisymmetric spacetimes. This result ensures 
us the equivalence in practice between the estimation of 
radiation exterior to the horizon and the verification of 
Penrose inequality via Dain's number. 



D. Multipolar analysis 

To be much more complete about the geometry of the 
constructed horizons, one could rely on the source multi- 
pole decomposition of the two-surface lying on our three- 
slice. This feature has first been presented by [65], based 
on an analogy with electromagnetism, and first studied 
in [37] in the case of dynamical horizons. We here im- 
plement the computation of multipole moments in the 
isolated horizon case, which is the strict situation where 
they have been defined in [65]. 

A prerequisite is the existence of a preferred 
divergence-free vector field (p a on the sphere, from which 
the angular momentum of the horizon is defined (the 
divergence-free condition on tp a ensures that all defini- 
tions will be gauge-independent). As mentioned above, 
our chosen vector field will be the one associated with 

the azimuthal coordinate, namely (^^j ■ 

Another important feature is the construction of a pre- 
ferred coordinate system, so that the Legendre polyno- 
mials associated with spherical harmonics will possess 
the right orthonormality properties; as expressed in the 
implementation of [37], this reduces to finding a set of 
coordinates (£, tp) where the metric on the two-surface 
can be written as: 



where A is the minimal area of a surface containing the 
horizon, Jk is the Komar angular momentum at infinity 



Qab 



Ru (/(c)- 



1 D a (D b C + f{£)D a <pD b <p) 



(38) 
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FIG. 4: Value of 1 — e.4 for both data sets. 

with R-j-c the areal radius of the sphere and /(£) deter- 
mined in terms of the two-dimensional Ricci scalar and 
the norm of ip a [65]. In the axisymmetric case studied 
here for the horizon, the integral curves for the coordi- 
nate ip are already defined by the orbits of the vector field 
(jp) 1 - The coordinate C is defined by 

= -i-ebaP 6 - (39) 
H n 

An appropriate normalization should be added, that 
ensures that § H C,d 2 V = 0. In the Kerr case, those coor- 
dinates turn out to correspond with the Boyer-Lindquist 
coordinates, with £ = cos 9 in spherical coordinates [37]. 

The mass and angular momentum multipolcs of order 
n are then defined, by analogy with electromagnetism 
[65]: 

cm i\ , r p 

Mn = j {nPn (Q}#V, (40) 

pn — l p 

J n = ^ j s P' n {QK ab s a V b d 2 V. (41) 

With this definition and using the Gauss-Bonnet the- 
orem it is trivial to see that Mq = Mh and J\ = Jk, the 
Komar angular momentum on the horizon. 

We should emphasize that these multipoles, except for 
Mq and J\, are in general different from the field gravi- 
tational multipoles that can be defined at infinity. How- 
ever, the authors in [65] have pointed out that the knowl- 
edge of all the multipoles of an isolated horizon allows to 
reconstruct the whole horizon, and also the spacetime in 
a vicinity of this horizon. The multipolcs then discrim- 
inate exactly every isolated horizon, and the spacetime 
at its vicinity. Figure 5 shows the capacity of telling 
apart the horizon of a CF axisymmetric slice and the one 
of a NCF slice, in Dirac gauge. Data are also compared 
with an analytic Kerr solution in Kerr-Schild coordinates. 
Apart from the accuracy obtained for our NCF data (and 
a further confirmation that we indeed have obtained the 
actual Kerr spacetime), we see the clear distinction made 
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FIG. 5: Computations of the rescaled second-order mass mul- 
tipole (Iff) and the rescaled third-order angular multipole 
(jp:)- Relative differences with respect to values for a Kerr- 
Schild analytical horizon are displayed. 

by this computation between the Kerr horizon and a con- 
formal approximation of it. Together with Dain's num- 
ber, this study has proven that those two tools are very 
well-suited to study isolated horizon properties, and the 
distance between data obtained from, say an evolution 
scheme, and the eventual equilibrium black hole data it is 
supposed to reach. Ultimate tests on the characterization 
of the obtained data as slices of Kerr could be achieved 
by implementing the schemes proposed in [66, 67]. 

VI. DISCUSSION 

The data we get with our simulations are interesting 
at several levels. They allow to make a direct compari- 
son between the conformally flat approximation and the 
exact solution for axisymmetric spacctimcs containing a 
black hole, that arc both calculated a priori and in the 
same gauge. As we have seen in Sec. V, this gives us in- 
sight about the geometric features of the exact solution; 
we can single important issues concerning for example the 
intrinsic geometry of the horizon, via multipoles and the 
Penrose inequalities. Numerical tools are in this respect 
implemented and their efficiency tested. 

At a more theoretical level, the method we used to 
get those data is a little bit heterodox: providing stan- 
dard non-expanding horizon conditions for (3+1) vari- 
ables such as j3 l , ip and Nip, we choose in addition not 
to prescribe any further geometrical information for the 
conformal part, symbolized here by the tensorial field h tJ . 
This has been motivated by the fact that, given the tenso- 
rial equation corresponding to h lJ in our formalism, it ap- 
pears that we most likely cannot prescribe anything else 
in the studied setting. By numerical transformation of 
the operator acting on h lJ we ensured that at every iter- 
ation step no boundary condition was required. The fact 
that our system of Einstein equations written this way 
converges to the required Kerr solution shows that in- 
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deed, no additional information was needed in the single 
horizon case. To be more precise, our method suggests 
that in this case, the conformal geometry of the MOTS is 
directly encoded in the Einstein equations, when written 
with an equilibrium ansatz: we directly use these equa- 
tions to justify a no-boundary treatment. 

In the light of this numerical study, we can make a par- 
allel with the proposition made in [39] . In that paper, the 
authors suggested, after a gauge dependence analysis for 
7 U on the horizon, that a prescription could be made for 
the conformal geometry on the horizon. In this respect, 
they justify, for the projection of the three-metric on the 
two-surface, the following choice: 

q ab = to 2 f abl (42) 

with f a b the usual diagonal round metric for a two-sphere 
in spherical coordinates adapted to the horizon. This 
choice should not affect the physics of the three-slice, 
and suffices to recover the solution for Einstein dynam- 
ical equation with a slice of a spacetime containing an 
isolated horizon. Their study is made in a differential 
gauge generalizing the Dirac gauge we are using, namely 
D il^ = V i with V 1 a regular vector field on the three- 
slice. Our case corresponds to V 1 = 0, which is precisely 
the one treated in detail in [39]. 

When comparing to our data, we find that the projec- 
tion of our three-metric on the two-surface is not con- 
formally related to the flat metric in adapted spherical 
coordinates. This means that in our particular case, and 
in regard of the particular no-boundary argument we use, 
a boundary condition of the same type as (42) is probably 
inconsistent with our data, and is likely a choice that we 
do not have the freedom to make (note however that ge- 
ometric conditions in [68] making full use of the isolated 
horizon structure are indeed compatible with the present 
results, i.e. they are identically satisfied in the present 
Kerr case, whose horizon is indeed an isolated horizon). 
Unfortunately, the authors in [39] did not present any 
numerical results to support their claims, that we could 
have compared with ours. 

We insist here on an important caveat for our argu- 
mentation: assumptions can only be justified in the very 
particular case we are studying here, which is the axisym- 
metric vacuum spacetime. This spacetime has very spe- 
cific and non-trivial properties, all related to the unique- 
ness theorem of Carter [69]. Although the reasoning we 
have made on Sec. Ill for the operator could apply in 
other isolated horizon studies, we are not certain that 
our algorithm would globally converge when applied to 
a more general case (e.g. a black-hole binary system); a 
failure of this behavior would probably mean that an ad- 
ditional information about the conformal two-geometry 
has to be given to the system. Geometric fully isolated 
horizon boundary conditions proposed in [68] could then 
be enforced (note that geometric inner boundary con- 
ditions in [68] are not necessarily tied to the particular 
analytic setting here discussed and, more generally, they 
would also apply in schemes not enforcing the coordi- 



nate adaptation to stationarity at the horizon, b = N-n, 
crucial for the singular nature of operators (30)). 

Finally, let us point out the fact that we made those 
simulations by prescribing only the geometry of the hori- 
zon, and the geometry of spacetime at infinity. No as- 
sumption has been made for axisymmetry in the three- 
slice (computations can easily be made in the full 3D 
case and give the same results). Prescribing a vanishing 
expansion and a conformal Killing symmetry on a hori- 
zon, together with asymptotically flat hypothesis, our 
code converges to the only solution of the Kerr space- 
time. Without claiming any rigorous demonstration here, 
this numerical result is most likely a support to the well 
known black hole rigidity theorem [28], where the same 
hypotheses lead to a uniqueness theorem involving the 
Kerr solution as the only one with no electromagnetic 
field. 



VII. CONCLUSION 

We have used the prescription for a fully constrained 
scheme of (3+1) Einstein equations in generalized Dirac 
gauge [25, 26] to retrieve stationary axisymmetric black 
hole spacetime, and compared it with the analytical so- 
lution of Kerr type. An advanced handling of the con- 
formal geometry of our three-slice allowed us to reach 
actual stationarity with good resolution precision for our 
scheme. Although we used standard quasi-equilibrium 
conditions concerning boundary values for other metric 
fields in the excised horizon, we found that the conformal 
geometry on the horizon required no prescription what- 
soever in the single horizon case. This is in contrast with 
suggestions available in the literature [39], and probably 
suggests an underlying physical feature of the horizon 
geometry (maybe related to uniqueness of the Kerr solu- 
tion). To our knowledge, it is the first time the conformal 
part is numerically computed in a black hole spacetime 
using only a prescription on the stationarity of spacetime 
(and without resorting to additional symmetries). The 
application of this feature to the more general initial data 
problem is evident: in the same spirit as the work done in 
[49] for neutron-star binaries, using it for the black-hole 
binary system could lead to significant improvement in 
the available initial data for evolution codes. Further 
numerical work will clarify this issue. 

We have implemented and used in our study numerical 
tools aimed at characterizing the geometry and physical 
properties related to horizons embedded in spacetime; 
those tools, among which a complete multipole analy- 
sis for two-surfaces as gravitational sources, have proven 
very accurate for diagnostics involving the horizon geom- 
etry and physical features. They will be more thoroughly 
presented, and tested in more general cases, in an upcom- 
ing work. 



13 



Acknowledgments 

The authors warmly thank Eric Gourgoulhon and 
Philippe Grandclcmcnt for many valuable discussions, 
and a careful reading of the manuscript. We also thank 
the referee for numerous useful comments. NV and JN 
were supported by the A.N.R. Grant BLAN07-1_201699 
entitled "LISA Science" . JLJ has been supported 
by the Spanish MICINN under project FIS2008-06078- 
C03-01/FIS and Junta de Andaluci'a under projects 
FQM2288, FQM219. 



APPENDIX A: TENSOR SPECTRAL 
QUANTITIES ADAPTED TO THE DIRAC 
GAUGE 

We here give the definition of the three spectral quanti- 
ties introduced in Sec. Ill C, that describe the divergence- 
free degrees of freedom (with respect to the Dirac gauge) 
associated with a rank two symmetric tensor. The reader 
is also invited to go to [52] or [53] where more detailed 
calculations are provided. 

We first define a set of spin-weighted tensor spheri- 
cal harmonics components for a symmetric rank-2 ten- 
sor, directly linked to the tensor spherical harmonics as 
introduced by Mathews and Zerilli [70, 71]. We shall 
give the expression for these components of the tensor 
h v using the classical spherical coordinate basis, which 
is used in practice in our computations. With the no- 
tation P = h ee + h vv , the six pure spherical harmonics 
components of h %1 are defined as : 
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d 2 h e v 



sin 2 e dp 2 
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the fifth and sixth scalar fields being simply the tensor 
trace h with respect to the flat metric and the h rr spheri- 
cal component. Let us note that these relations are more 
tractable when using a scalar spherical harmonics decom- 
position (introduced in Sec. IV B) for all fields. Indeed, 
an angular Laplace operator acting on a field reduces 
then to a simple algebraic operation on every spherical 
harmonics component. Inverse relations can also be com- 
puted to retrieve the classical components of h %3 from 
spherical harmonics quantities. 



We now derive the main variables related to our study: 
with the divergence- free decomposition h %3 = T> l W 3 + 
T> 3 W l + hj., and T>ihj. = 0, a choice for three quantities 
defined from h lJ and verifying: 



lit 3 = 



A = B = C = 0, (A5) 
can be expressed as the following scalar fields (see [53]): 
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These quantities can also be decomposed onto a scalar 
spherical harmonics basis. The equivalence in (A5) is 
achieved up to boundary conditions. 

To show how the quantities A, B and C behave with 
respect to the Laplace operator, we shall assume in the 
following that the tensor h %3 is the solution of a Poisson 
equation of the type Ah 13 — S 13 . We can deduce a scalar 
elliptic system verified by A, B and C as: 



AA = A s 

AB -^ = Bs 
AC + ^ + ^ = C s , 



(A9) 
(A10) 

(All) 



Where As, B$ and Cs are the corresponding quantities 
associated with the source S 13 . A simple way of decou- 
pling the last two elliptic equations is to define the vari- 
ables B = Elm B lm Yim and C = E ; m & m Y lm with: 



B' 
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Thus, we can write an equivalent system 
for (A9,A10,A11) as: 



AA = A s , 
AB = B s , 
A*C = C s , 



(A14) 
(A15) 
(A16) 



With the following elliptic operators defined for each 
spherical harmonic index 1: 
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/ is the identity operator. A, B and C arc then defined 
as three scalar fields characterizing only the divergence- 
free part of a symmetric rank two tensor h %3 , and giving 
a system of three decoupled scalar elliptic equations in 
the Poisson problem for this tensor. Hence they are very 
well suited to the study of a tensorial elliptic problem in 
Dirac gauge. Let us finally note that the quantities B 
and C are directly related to each other by the trace of 
the considered tensor [53] : If we know a priori the value 
of the trace for h 13 , then the knowledge of B suffices to 
recover C with no additional information (the converse 
being equally true). 



APPENDIX B: RECOVERY OF ti 1 FROM A AND 

B 



the quantities A and B, will allow to express two de- 
coupled differential systems. The first one, involving the 
spherical harmonics components fi and X , combines the 
expression for the scalar field A, as well as the fact that 
vanishes under the Dirac gauge : 



dX 
~dr~ 

du in, 1 , A 
or r r 



(B5) 



2)X = 0. 



The second system is composed of the definition of B 
for each of its spherical harmonic component B , as 
well as the vanishing of H r and H v (again, due to Dirac 
gauge): 



In this section, we come up with technical details for 
resolution of the gauge differential system introduced in 
Sec. Ill C, to reconstruct the tensor h %3 from the quanti- 
ties A and B. 

We begin by expressing components of the vector field 
representing the divergence of h 13 : 



IT = 'P, //'■'. 



(Bl) 





= B lm + 


dh rr 


3h rr 
H h 


dr 


r 


dn 


32 + l( 


7T 4 
or 


r r \ 



r 



2(1 + iy 

- (Ag^r] -h) = 0, 
r 



h - h rr 



(B6) 



= 0. 



with the Dirac gauge for H % = 0. Adopting the vec- 
tor spherical harmonics decomposition suggested in [25], 
the three spherical harmonics components of H* are ex- 
pressed, in function of the spherical harmonics compo- 
nents of h tJ (see Appendix A), as: 



W 



dh r 



3h rr 1 



in = a 9v 

H* = Ag^ 



(B2) 



^ + ^+ 1 -(A 6v + 2)X 
or r r 



(B4) 



Those three expressions, alongside with definitions of 



with the expressions (A7, A8) of B and C as functions of 
the spherical harmonics components of h lJ . In this sys- 
tem, the trace is given a priori, so that only three spher- 
ical harmonics components are considered as unknowns. 

We refer to the analysis of [53] to affirm that, when 
solving our equations in M 3 minus an excised inner 
sphere, one boundary condition has to be provided at 
the surface for the system (B5), and two for the system 
(B6). As pointed out in Sec. IV B, these conditions are 
retrieved as compatibility conditions based on the orig- 
inal elliptic tensorial equation. Overall, we are able to 
invert the two Dirac differential systems, and retrieve all 
the spherical harmonics components of h tJ from the sole 
knowledge of A, B and the trace h. 
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